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We present a model of adaptive regulatory networks consisting of a simple biologically-motivated 
rewiring procedure coupled to an elementary stability criterion. The resulting networks exhibit a 
characteristic stationary heavy-tailed degree distribution, show complex structural microdynamics 
and self-organize to a dynamically critical state. We show analytically that the observed criticality 
results from the formation and breaking of transient feedback loops during the adaptive process. 
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Introduction - Much recent research attention has fo- 
cused on understanding the structure of naturally oc- 
curring empirical networks and associated random graph 
models [31 0] . An overarching aim of many of these stud- 
ies is to determine the relationships between network 
structure and dynamics. For instance, the presence of 
modularity and sparsity have long been known to con- 
tribute to global stability, while the presence of feedback 
is a well-studied prerequisite for the support of complex 
dynamics such as oscillations, multistability and chaos 
[3(6]. However, although much work has, so far, focused 
on networks which are static in their topology, many real- 
world complex systems evolve both structurally and dy- 
namically over time [7HH]- For instance, neural networks 
change in structure depending on synaptic activity while 
genetic regulatory networks change structurally on the 
evolutionary time scale in a fitness-dependent manner. 
Consequently, adaptive networks - in which changes in 
network topology and dynamics continually feedback on 
each other - are now attracting increasing research in- 
terest jTUj. Since many biological regulatory systems, 
such as neural and genetic regulatory networks, are also 
thought to optimally balance stability and adaptability 
by operating at, or near to, criticality |11H14| . a num- 
ber of prior studies have sought to elucidate mechanisms 
of self-organized criticality (SOC) in adaptive networks 
[15H22] . For example, important early results were ob- 
tained by Christensen et al. and Bornholdt and Rohlf, 
who showed that adaptive networks may self-organize to 
a critical state by a simple mechanism in which 'quiet 
nodes grow links [and] active nodes lose links' [231 124j. 
However, despite the apparent ubiquity of critical adap- 
tive networks in nature, the mechanisms of adaptive SOC 
remain to be fully determined. 

In this article we outline a simple new adaptive net- 
work model which reproduces characteristic features of 
biological systems, including a heavy-tailed degree dis- 
tribution and self-organization to a dynamically critical 
state. To fix ideas our model may be thought of as 
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describing adaptive changes in a genetic regulatory net- 
work, although the model may also be applied more gen- 
erally to other systems which undergo adaptive rewiring. 
In genetic regulatory networks, genetic mutations cause 
changes in protein structure which, in turn, not only alter 
local network connectivity but also global system stabil- 
ity. Consequently, our model is a simple scheme intended 
to describe, albeit in a highly idealized way, mutation- 
driven local rewiring in the face of a global stability (fit- 
ness) constraint: mutations are allowed to accumulate 
during times of stability, but harmful mutations are sup- 
pressed during times of instability. 

Preliminaries.- Mathematically a network is a graph 
consisting of a set of vertices (or nodes) V of size n and a 
set of edges (or links) E. A directed graph (digraph) is a 
graph in which each edge Vi ~ vj £ E has a unique orien- 
tation (vi Vj). Although digraphs describe well struc- 
tural relationships in complex systems, in many cases 
relationships also have an intrinsic sign - friendship and 
enmity in social networks or activation and inhibition in 
biochemical regulatory networks, for instance. To cope 
with such systems, a natural framework is that of signed 
digraphs. A signed digraph S is a digraph in which 
each edge Vi ~ vj 6 E additionally has a unique sign 
<Tij G [— 1, +1] depending on whether it is 'activating' 
(<7jj = +1) or 'inhibiting' (try = —1). The adjacency ma- 
trix A = a+j of a signed digraph has the form ay = cry 
if Vi ~ vj € E and ay = otherwise. When considering 
structural features of S without regard for signs we shall 
also make use of the absolute adjacency matrix A = |ay |. 
The in-degree (out-degree) of a vertex is the number of 
in-coming (out-going) edges it has, without regard for 
sign. The net-degree d net (vi) = \d in (vi) - d ou t(vi)\ of a 
vertex Vi as the absolute difference of its in-coming and 
out-going degree. Intuitively, net-degree measures how 
'source-' or 'sink'-like a vertex is. By extension, we define 
the imbalance of a vertex-pair as the absolute difference 
of their net-degrees, I(vi,v v ) = \d DBt (vi) - d net (Vj)\. It 
has recently been observed that many empirical networks 
contain significantly more source- and sink-vertices than 
expected by chance, and that this degree imbalance nat- 
urally leads to depletion of feedback loops (cycles) which, 
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Figure 1: Graphs resulting from the evolutionary process ex- 
hibit a stationary heavy-tailed degree distribution and com- 
plex microdynamics. (Left) The net-degree distribution at 
1000 time-step intervals for a period of 1 x 10 6 time-steps at 
equilibrium for a typical 250 vertex network are shown in light 
gray. In bold (red online) is the mean net-degree distribution 
over this time period. (Right) Each row shows the changing 
net-degree of one vertex over the same time period. 



in turn, confers enhanced stability properties |25j . Thus, 
degree imbalance and dynamic stability are intrinsically 
related, a fact that our model exploits. 

Model- We begin at t = with a random signed di- 
graph S(t = 0) of size n with Erdos-Renyi connectivity, 
in which edge orientations and signs have been assigned 
independently in an equiprobable random manner |41j . 
We then rewire S(t) at successive time-steps according 
to the following rules: (1) randomly and uniformly chose 
an edge e id = v a ~ Vb connecting two vertices in S(t) 
such that S(t) — e id is not disconnected and an ordered 
pair of non-adjacent vertices v c ,Vd ^ v c . (2) Calculate 
the pair-wise imbalances l(v a , Vf,) and I(v c , V4). (3) Delete 
e id and create a new edge e ncw = v c — > va, choosing its 
sign randomly and uniformly, and recalculate the imbal- 
ances. (4) If the sum of the two imbalances after the 
switch is greater than that before then accept the switch 
unconditionally, otherwise accept with probability p(t). 

In order to couple structural rearrangement to dynam- 
ics we allow p(t) to vary in a manner which takes into 
account the changing stability of the system. To do so 
we assume that, in addition to regulatory links defined 
by S(t), each species (vertex Vi) also decays at a constant 
characteristic rate which we fix at t — independently, 
randomly and uniformly on the unit interval. Thus, at 
each evolutionary time-point we obtain a modified adja- 
cency matrix B(i) = A(t) — diag(ej). Global stability is 
then given by the magnitude of /z max (£) = max Re Pi(t), 
where Hi(t) for i = l...n are the eigenvalues of B(i). 
In particular, the system is stable when /i max (i) < 
and unstable when /i max (£) > [42 . Therefore we set 
p(t) = 1 — h{p maK (t)], where h[x] is the Heaviside step 
function, allowing defective switches when the system is 
stable and suppressing defective switches when the sys- 
tem is unstable. The key property of this coupling is that 
it makes global information available to the local struc- 
tural reorganizing process, providing continual feedback 
between structure and dynamics. 

Results - The networks produced by this simple model 



are characterized by a stationary heavy-tailed degree dis- 
tribution (see Fig. [Tjleft) indicating the presence of hub 
source- and sink- vertices, a well-known feature of real- 
world networks [351 HI]. However, since our model al- 
lows for periods of random structural rearrangement, this 
macroscopic stationarity masks complex structural mi- 
crodynamics in which individual vertices continually ac- 
cumulate and lose edges and rise and fall in their cen- 
trality (see Fig. [I] right). This kind of 'mixing' micrody- 
namics is not produced by classical rich-get-richer mod- 
els of hub formation [25], but has recently been high- 
lighted as an important characteristic of real-world evolv- 
ing (macroscopically stationary) complex networks [7HS] . 

Fig. [2] gives a plot of /x max (£) at equilibrium [33] for a 
representative system showing that dynamics on the evo- 
lutionary time-scale are characterized by periods of sta- 
bility (/x max (i) < 0) and instability (/U max (i) > 0) punctu- 
ating back- and- forth. To help interpret these dynamics, 
also shown is A max (t) = max Re Ai(i), where Xi(t) are the 
eigenvalues of the graph adjacency matrix A(i) and three 
measures of network structure. The first structural mea- 
sure shown is total net-degree D nct (t) = d net [vi(t)], 
a measure of overall degree imbalance in S(t). It is ap- 
parent that changes in total net-degree correlate poorly 
with changes in stability, suggesting that although fluc- 
tuations in net-degree are observed during the evolution- 
ary process, it is not degree-imbalance per se that drives 
the characteristic dynamics of £t max (i)- In order to iden- 
tify more precisely the structural origin of the observed 
bursting dynamics, and based upon the observation that 
degree imbalance naturally leads to feedback loop deple- 
tion [2 5) , also shown are two measures of network cyclic 
structure [H]. The first, Q(t) — n cyc (t)/n where n cyc (t) 
is the number of vertices which participate in a cycle 
in S(t), measures overall cyclic structure without regard 
for details such as cycle numbers or distribution of cycle 
lengths. The second, 

n 

(i) = Trace e A(t) - n = ^ e^ - n, (1) 

where A^ are the eigenvalues of the absolute adjacency 
matrix A(i), is an indirect measure of 'returnability' 
which takes into account details of closed walks in S(t). 
In particular, ^(t) is a sum of all closed walks in S(t) 
weighted in decreasing order by length (note that ^(t)+n 
may be thought of as the partition function of S(t)) 

[3HE2]. 

Examining the time-series of vf'(t) and $(£) it is ap- 
parent that, unlike total net-degree, both ^f(t) and $(t) 
exhibit similar bursting behavior to that of /x max (t). In 
particular, periods of stability (/^ max (£) < 0) generally 
correspond to periods when both *ff(t) — and $(<) = 
(in Fig. [2] this occurs > 90% of the time). Since <&(i) = 
and (f>(t) = if and only if S(t) is acyclic this indicates 
that periods of stability occur primarily when S(t) is 
acyclic. Furthermore, changes in stability predominantly 
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Figure 2: Evolutionary dynamics are characterized by punc- 
tuated equilibrium. Dynamics of a representative 100 vertex 
network are given. (Top) shows the time-series for /i max (t). 
Note that /i ma x(t) = — mine 4 < when S(t) is acyclic (in 
this case mine^ = 6 x 10~ 3 ); (Second) the time-series for 
Amax(t), note the prevalence of and +1 (and less obviously 
to the eye, but still present, cos(7r/Z) for some I € Z + ) in this 
series as predicted analytically; (Third) the total net-degree 
D nc t(t) = J2idnet[vi(t)]; (Fourth) the cyclic index <&(£); (Bot- 
tom) the cyclic index \P(t). 



occur concordantly with changes in \I>(t) and $(i) (for 
instance, in Fig. [2] this occurs > 99% of the time). Con- 
sidering the time-series as binary variables ('stable or un- 
stable' and 'cyclic or acyclic') and calculating entropies 
gives #[$(*)] = H[V(t)] = 0.093 and ff^ft)] = 
0.100. The mutual information between these series is 
M[/w(*),$(*)] = MDwW,*®] = 0.088, indicating 
that changes in stability are strongly, although not ex- 
clusively, related to changes in cyclic structure (see also 
Fig.§ 05]. 

In order to better understand this relationship, we now 
derive some analytical results relating cycles and spectra 
of signed digraphs which will help interpret these numer- 
ics. To obtain exact results we shall focus on deriving 
analytic formulae for "J and A max in the particular case 
that all cycles in S are disjoint (that is, each vertex v E V 
belongs to at most one cycle). Although this is a strong 
condition to impose, and most real- world networks are 
not expected to be cycle-disjoint, this case is analytically 
tractable and, since our evolutionary scheme favors the 
minimization of cycles, yields results which shed light on 
the observed dynamics. Full proofs of all analytic results 
are provided in Appendix |A"| 

Firstly we observe that if a signed digraph S is cycle- 
disjoint, then its spectrum has a particularly simple form. 
Specifically, if S contains c£ positive cycles and c^T nega- 
tive cycles of length fc (for k = 3 . . . n) [H] and all cycles 
are disjoint, then its spectrum is the zero eigenvalue with 
multiplicity (n~ n cyc ), along with the eigenvalues of each 
of the cycles considered separately as induced subgraphs 
(that is, the union of c£ copies of the fc-th roots of +1, 
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Figure 3: Transient cycles trigger bursts of instability. (Left) 
A plot of A*(i) = 9(t) - 9(t - 1) against A/i max (t) = 
fi m ax(i) — Mmax(i — 1) using the same data as Fig. [2] For 
clarity, the few changes in stability which do not occur con- 
cordantly with changes in cyclic structure are shown in light 
gray (red online). (Right) A close-up of the left panel. The 
striations arise since it is common for isolated cycles to be 
created or broken during the evolutionary process, triggering 
changes in stability. The dotted vertical lines are at ±^ cal- 
culated analytically using Eq. [2] with k = 3 . . . 6 and Ck = 1. 



and c fc copies of the fc-th roots of — 1, for fc = 3 . . . n). 

An immediate consequence of this result is that if S 
is cycle-disjoint and possesses at least one positive cy- 
cle then A max = 1, while if all cycles are negative then 
Amax = Re e 7 ™/' = cos(tt/1), where I is the length of 
the longest cycle in S. In this sense positive cycles are 
uniformly destabilizing, while the destabilizing effect of 
negative cycles increases with length. Examination of the 
time-series data shows that A roax (i) = 0, 1 and cos(n/l) 
for some 3<I<n6 Z + , do indeed occur commonly dur- 
ing evolution (for instance, in Fig. 2] this occurs ps 53% of 
the time), indicating the continual formation and break- 
ing of isolated cycles by the evolutionary scheme. 

This result is also useful since it allows us to calculate 
\E' analytically in the case that S is cycle-disjoint. In 



particular, if S contains Cfc (= 
length fc for k = 3 . . . n 



E 

fc =3 



then, using Eq. [T] 
Cfc H k ,o(l) - n c 



') disjoint cycles of 



(2) 



where Hk,o(z) is the generalized hyperbolic function of 
order fc and kind 36J. Fig. [3] shows that values of ^(t) 
calculated using Eq.[2]often occur during evolution, again 
indicating that isolated cycles are continually formed and 
broken by the evolutionary scheme. 

These analytical results may be used to interpret nu- 
merics by making use of two further results which relate 
Amax to /i max in the cycle-disjoint case. Firstly, note that 
in the special case that S(t) is acyclic then A max (i) = 
and /i max (*) = — minej < 0, and the system is stable. 
Secondly, if S(t) is cycle-disjoint then ^ max (t) < A max (t) 
and this bound is tight (/x max (t) — > A max (i) as — > 
for all i). Consequently, if vertex decay rates are all 
small then the completion of a single cycle in an oth- 
erwise acyclic network is sufficient to trigger a burst of 
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Figure 4: Bursts of instability have heavy-tailed statistics. 
The distribution of burst durations is shown over an interval 
of 4 x 10 6 evolutionary time-steps for a representative system 
with 100 vertices at equilibrium. A power-law with exponent 
1 is also shown for reference. 



instability, as seen in Fig. [2] When this occurs the evo- 
lutionary process responds by suppressing any further 
defective switches and rearranging local network struc- 
ture to remove the cause of the instability. Typically, 
this is quickly achieved and the burst of instability is 
relatively short. However, occasionally cycles may accu- 
mulate more rapidly than they are removed, giving rise to 
extended bursts of instability and heavy-tailed statistics 
characteristic of a critical state (see Fig. [4]) . 

For completeness it should be noted that if S(t) is not 
cycle-disjoint then the relationship between cycles and 
stability can be complex: it is not necessarily true that 
/■tmax(i) < A max (i) and, in rare cases, changes in cyclic 
structure and stability may occur discordantly. Further 
details of when this occurs and a simple worked-example 
are included in Appendix [B| 

Conclusions.- Many biological regulatory systems are 
thought to balance stability and adaptability by self- 
organizing to a dynamically critical state [TTMl4| . In 
this article we have presented a simple adaptive net- 
work model which reproduces characteristic features of 
biological systems, including a heavy-tailed connectivity 
distribution, microdynamics and robust self-organization 
to criticality. Previous models have shown that adap- 
tive networks may self-organize to a critical state due 
to rewiring based upon local activity [T5T424j . Here, the 
mechanism of self-organization is somewhat different and 
relies on the fact that feedback and stability are gener- 
ally inversely related: by employing a flexible rewiring 
scheme which allows feedback loops to be formed dur- 
ing periods of stability and eliminated during periods of 
instability, criticality naturally arises in our model. It 
seems plausible that these (and other, as yet unknown) 
adaptive processes may be responsible for the criticality 
observed in nature. 



Appendix A: Proofs 

Formal statements and proofs of results stated without 
proof in the main text are provided here. 

Let S = S(V,E) be a directed signed graph with 
vertex-set V of size n, edge-set E of size m, no dual- 
edges (if ViVj e E then VjVi ^ E for alii, j = 1 . . . n) and 
no self- loops (viVi ^ E for alH = 1 . . . n). We say that 
S is cycle-disjoint if all cycles in S are pair-wise disjoint, 
that is, if every vertex participates in at most one cycle. 
Let c^T (respectively c^T ) be the number of positive (re- 
spectively negative) cycles of length k in S for k = 3 . . . n. 
Thus, the number of vertices which participate in a cycle 
in S is n cyc = J2k k (4 + c k)- 

Proposition 1. The eigenvalue spectrum of a cycle- 
disjoint directed signed graph S consists of the zero eigen- 
value with multiplicity (n — n cyc ) along with the union of 
c^ copies of the k-th roots of unity and c^ copies the k-th 
roots of —1 for k = 3 . . . n. 

Proof. The proof makes use of Sachs' (coefficients) theo- 
rem (Theorem 1.32, p32 in [37]) which, for completeness, 
we state here in its general form. 

Theorem 1 (Sachs). Let W be a weighted digraph with 
characteristic polynomial z n + a\z n ~ l + . . . + a„_iz + a n 
then 

= E (-1) P{L) W(L), (Al) 

LeC z 

where Ci is the set of directed linear subgraphs of W on 
i vertices, p{L) is the number of disjoint components in 
a given linear subgraph L and W{L) is the product of 
edge-weights over all edges in L. 

We now begin our proof of Proposition [l] Let A be 
the adjacency matrix of a cycle-disjoint signed digraph 
S. The eigenvalues of S are the solutions to the charac- 
teristic polynomial of A 

z n + aiz"" 1 + . . . + a n -iz + a n . (A2) 

The largest linear subgraph L max in S consists of the 
disjoint union of all the cycles in S and so has size n cyc = 
Sfc k( c t + c fe )■ Therefore, by Sachs' theorem ai — for 
i > n cyc and the eigenvalues of S are solutions to 

z n + a x z n - x + . . . + a ncyc z n ~ n ° yc = (A3) 
+ aiz"^- 1 + . . . + a ncy J = 0. (A4) 

Thus, S has zero as an eigenvalue with multiplicity (n — 
n cyc ). Now let 

z n ^ + a x z n <v°- 1 + ... + a ncyc _ lZ + 5„ cyc (A5) 

be the characteristic polynomial of L max , considered as 
an induced subgraph. It is immediate from the definition 
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of L max (and using Sachs' theorem) that 5j = a, for all 
i. Thus, the additional eigenvalues of S are the roots of 



Eq. 



5] which are the eigenvalues of the disjoint cycles in 
and the result follows. □ 



Proposition 2. Let S be a cycle- disjoint signed digraph 
with n vertices, adjacency matrix A and eigenvalues 
Ai, . . . , A„. Let b = diag(—ei, . . . , — e„) be a diagonal 
matrix with > for all i and let //!,...,//„ 6e t/ie 
eigenvalues of the matrix B = A + b. Then 

maxRe(^) < maxRe(Ai) 

i i 

with equality only when ti = for all i. 

Proof. The proof consists of three parts: (1) reduction of 
the problem to that of a cycle; (2) proof for the positive 
cycle case; and (3) proof for the negative cycle case. 

Part 1: reduction of the problem to that of a cycle. 
The matrix B = A + b may be thought of as the ad- 
jacency matrix of a weighted digraph P which has the 
same vertices, edges and edge-signs as 5* with an addi- 
tional self-loop at each vertex Vi £ P of weight — ej. Let 
v be a vertex not participating in a cycle in S. Any lin- 
ear subgraph L of P containing v can only do so via the 
self-loop at v and hence v must be disjoint from all other 
vertices in L. Thus, removing all edges in P which do not 
participate in a cycle in S, except the weighted self-loops, 
creates a new graph Q with the same linear subgraphs 
as P and thus the same characteristic polynomial as P 
by Sachs' theorem. The graph Q consists of the disjoint 
union of the perturbed cycles of S and (n — n cyc ) isolated 
vertices each with a self-loop of weight — for some i. 
Each of the isolated vertices contributes the eigenvalue 
— 6i < for some i to the spectrum of Q. The remainder 
of the spectrum is determined by the disjoint union of 
perturbed cycles in Q. Since the spectrum of a disjoint 
union of graphs is the union of the spectra of each of the 
components, it is therefore sufficient to prove the propo- 
sition for S a cycle. In particular, if S is acyclic then 
the spectra of P is just the set {— | i = 1 . . . n} which 
has maximal real-part — min q < as stated in the main 
text. 

Part 2: proof for the positive cycle case. Let G+ 
be a positive cycle with n vertices and adjacency ma- 
trix A. The characteristic polynomial of C+ is z 11 — 1 
and the eigenvalues of are therefore the nth roots 
of unity, which have maximal real-part 1 for any n. 
The adjacency matrix of a perturbed positive cycle is 
A + diag(— ei, . . . , — e„), which has characteristic polyno- 



mial (cf. Eq. Al I 



p+{z) = (z + e 1 )-(z + e 2 )-...-(z + e„)-l = JJ^+e^-l. 



Thus, we need to prove that every root A of p+ satisfies 
Re(A) < 1 when at least one ti £ M. is nonzero. From 



now on we shall assume, without loss of generality, that 
ei > 0. 

We shall use Rouche's theorem, a well-known theorem 
in complex analysis for locating the roots of functions. 
For a proof of Rouche's theorem see [35]. 

Theorem 2 (Rouche). Let f and g be holomorphic func- 
tions in a domain R C C. Let D C R be a bounded subset 
such that its boundary dD is a simple closed curve null- 
homologous inside R. If 

\f(z) - g(z)\ < \f(z)\ for all z E dD (A6) 

then f and g have the same number of zeros inside D. 
To apply Rouche's theorem we take 



9(z) 



— «+ 



Pn( z ) 



and we will define D in a moment. It is immediate that 
\f(z)~ g(z)\ — 1 for all z € C. Let be the ball of radius 
1 centered at — ej for each i and £8 = [j i be the union 
of the balls. Then 



!/(*)! = 11 + ifz ^ 



(A7) 



Now f(t) with t € [0, 1] is a strictly increasing real func- 
tion with /(l) = n™=i(l + e i) > l + > 1 so there exists 
< t < 1 with /(t ) > 1. Indeed 

> 1 for all z e C with Re(z) = t . (A8) 

To see this, let z G C with Re(z) = to and Im(z) = y. 
Then 

\z + ti\ = V(to + ^) 2 + y 2 > Vito + e;) 2 = \to + ti\ 

for alH = 1, .. . ,n and thus \f(z)\ > \f(t )\ > 1. 
Finally, let a > 1 and b < mini{— — 1} and define 

D = {zeC\b< Rc(z) < t and - a < Im(z) < a}. 

Then for all z € dD we have either Re(z) — to or z ^ , 



so \f{z)\ > 1 = \f(z) - g(z)\ by Eq. A7 and Eq. A8 



Therefore Rouche's Theorem gives that f(z) and g{z 
both have all their roots inside D. In particular, any 
root A of g(z) — p+(z) satisfies Re(A) < t < 1 and this 
proves the positive case. 

Remark 1. Note that, in addition, any root of Pn{z) 
with positive real-part must lie strictly inside the unit 
circle. To see this, observe that if \z\ = 1 then \f(z)\ = 
Yli \ z + e i\ > k + > 1 so z is not a root of p+(z) by the 
same argument. We shall make use of this observation in 
a moment. 

Remark 2. The positive case may also be proven by 
application of a modification of Gersgorin's disc theorem 
due to Brualdi (see Theorem 6.4.18 in [55]). 
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Part 3: proof for the negative cycle case. Let C n be a 
negative cycle with n vertices and adjacency matrix A. 
The characteristic polynomial of C~ is z n + 1 and the 
eigenvalues of C~ are therefore the nth roots of —1, which 
have maximal real-part cos(7r/n). The adjacency matrix 
of a perturbed negative cycle is A + diag(— e±, . . . , — e„), 



and the result is proven. To do so, we make use of the 
following two lemmas, whose proofs we leave to the end. 



which has characteristic polynomial (cf. Eq. Al ) 



p-(z) = {z + e 1 )-(z + e 2 )-...-(z + e n ) + l = + + 1. 

i=l 

Thus, we need to prove that every root A of p~ satisfies 
Re(A) < cos(7r/n), when at least one e, € K is nonzero. 
Again assume, without loss of generality, that t\ > 0. 

We first note that by exactly the same argument as the 
positive case, we can prove that all the roots A of p~ (z) 
satisfy Re(A) < 1 and (by Remark 1) that any root of 
p~ (z) with positive real-part must lie strictly inside the 
unit circle. However, in the negative case this bound is 
not sufficiently tight to prove the result since the magni- 
tude of the maximal real-part of the eigenvalues depends 
upon the length of the cycle. In fact, we require the tight- 
est possible bound and the proof in the negative case is 
correspondingly more involved than that of the positive 
case. 

We shall progress as before. However this time we use 
a strengthened version of Rouche's Theorem (3^1 p. 390] 



in which the inequality in Eq. A6 is replaced by the in- 
equality 

\f(z) - g{z)\ < \f(z)\ + \g(z)\ for all z G 3D. (A9) 
In this case, we use the functions 

f(z) = z n + l, 

g( z ) = Pn( z )- 

and the region 

Dg = {z = re lip €C|0<r< 1,-0 <</><#} 
where < 9 < ir/n. In particular, we shall prove that 

\f(z) - g(z)\ < \f(z)\ + \g(z)\ for all z € dDg 

from which it follows that Pn( z ) has the same number 
of roots that f(z) in Dg. Since, by construction, f(z) 
has no roots in Dg this implies that all roots A of p~ (z) 
with positive real-part must lie strictly in the unit circle 
excluding the region Dg for all < 9 < ir/n and, in 
particular, that each root A of p~(z) has real-part less 
than cos(7r/n). 

First note that in general 

\f(z)\ + \g(z)\ = \f(z)\ + \g(z)-f(z) + f(z)\ 
> \f{z)\ + \g(z)-f(z)\-\f(z)\ 
= W)~f{z)\ 

so the non-strict inequality holds for every z € C (observe 
that equality holds, for instance, for any root of either / 
or g). We therefore only need to demonstrate that 

\f(z) - g(z)\ ? \f(z)\ + \g(z)\ for all z G dDg (A10) 



Lemma 1. Let u,v G C. Then 

\u\ = \u + v\ + \v\ 



(All) 



if and only if either u — v = or v — a u with < a < 
-1. 

Consider the (open) upper and lower half-planes 

H+ = {zeC \ Im(z) > 0}, 
H~ = {zeC \ Im(z) < 0}, 



and the four (open) quadrants 



Q k = { z = re^ e C | r > 0. 



- 1W kn , 



for k — 1 ... 4. Additionally, for a nonzero complex num- 
ber w write Arg(u>) for the unique (p G (— 7r, n] such that 
w = |iu|e^. 

Lemma 2. Let beR+. 

1. If z £ Qx then \z+b\ > \z\ > andO < Arg(z+6) < 
Arg(z). 

2. Suppose that w±,W2 G H + satisfy \w\\ > \vj2\ and 
Arg(wi) < Arg(u>2). Then either w± -11)26 H + or 
W\ — u>2 G Q 4 UR + . 

We now apply Lemma[T]and Lemma[2]to complete the 
proof. 

Divide Dg into three regions: 

Dj = Dg n H + , Dg = Dg n H~ and Dg n E. 

Let z G D~g . To make use of Lemma[l]set u — s n _iz n ~ 1 + 
. . . + S1Z + S0 and v — z n + 1. We argue by contradiction. 
Suppose that z does not satisfy Eq. |AlO[ that is, in terms 
of u and v 

|tt| = |tt + «| + |t;|. (A12) 

Therefore by Lemma [TJw and v lie on a line through the 
origin. We shall prove that v G Qi and u g 1 Q3 and hence 
arrive at a contradiction. 

We know that for z G Dg, < Arg(z) < 9 < n/n and 
\z\ < 1 hence < Arg(z") < it and \z n \ < 1. Conse- 
quently v = z n + 1 € Qi. On the other hand, consider 



Sn-lZ 



+ ... + Sl z + s = Y[(z + hi) 



Write w\ — II"=i( z + ^) an d w 2 = z n . By Lemmara(l), 
\w\\ > \v)2\ and < Arg(wi) < Arg(w 2 ) since at least 
bi > 0. By Lemma [2] (2), u = Wi — W2 lies in either H + 
or Q4 U 1R + and thus u £ Q3. 

If z G Dg , we apply exactly the same argument to the 
complex conjugate z G Dt to conclude that u and v do 
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not lie in a line through the origin, therefore neither do 
u and v, again contradicting Lemma [l] 

Finally, if z £ Dg ni = [0,1] then u and v are both 
positive real and hence do not satisfy Lemma [T] and this 
completes the proof. □ 

Proof of Lemma^ Write u = a + bi, v = c + di. Then 
|u + v\ = \v\ — \u\ means 

y/(a + c) 2 + (b + d) 2 = yfc? + d 2 - yfa 2 + b 2 

which implies that (squaring and simplifying) 

ac + bd= v/(a 2 + b 2 )(c 2 + d 2 ). 

Squaring and simplifying again we obtain 

(ad — be) 2 — 0, that is, ad = be . 

If it = Q then Eq. |Al~l~1 impl ies v = . If u 7^ then either 
a = d/b or a = c/a is well-defined and satisfies v = au. 
In addition, 

\u\ = \u + v\ + \v\ = \u + au\ + \au\ =£• l = |l + a| + |a| 

and a case study shows that a < and 1 + a > 0, that 
is, — 1 < a < 0. One finally checks that, for such an a, 
v = au satisfies Eq. \KTT\ □ 

Proof of Lemma^ (1) Let z = x + iy with x, y > 0. 
Then 

\z + b\ = y/(x + b) 2 + y 2 > yjx 2 + y 2 = \z\ 

since x + b > x > 0. Recall that arctan is a strictly 
increasing function. Thus 



Arg(z + b) = arctan 




(2) Let ttfi = xi + iyi and w 2 = x 2 + iy 2 - If Vi > Vi then 
\miw\ — w 2 ) — yi — y 2 > so w\ — w 2 £ H + . Suppose 
that y 2 > y\ > (the reader should draw a picture at 
this stage to convince themselves). Since \wi\ > \w 2 \ 
we have x\ — x 2 > y 2 — yf > 0, that is, x\ > x\ or, 
equivalently, \x\\ > \x 2 \. Then either (a) x\ > \x 2 \ or 
(b) x\ < —\x 2 \. The latter is impossible: if x 2 > then 
< Arg(w 2 ) < it/2 but x\ < -x 2 < so Arg(wi) > tt/2; 
if x 2 < then x\ < x 2 < and hence 

1 1 Vi Vl , 
— < — — < — =4> Arg(w 2 ) < Arg(wi). 
x 2 X\ x 2 x\ 

In either case the condition Arg(iui) < Aig(w 2 ) is contra- 
dicted. So we must have (a) x\ > \x 2 \, that is, X\— x 2 > 
and therefore 

Wi — w 2 = (x\ - x 2 ) + (yi - y 2 )i £ Q 4 U R + . □ 




Figure 5: A graph which exhibits dissipation-induced insta- 
bility. Ghost edges show how this graph may be embedded in 
a larger graph. 



Appendix B: Dissipation-induced instability 



Occasionally in our model positive and negative cycles 
will intersect in a locally symmetric manner such that 
their contributions to the graph spectrum cancel each 
other out. In these cases, disparate decay rates may act 
to break the symmetry, giving rise to 'dissipation-induced 
instabilities' [30] • A simple example of when this occurs, 
consisting of a positive and a negative cycle arranged 
back-to-back, is given in Fig. [5] In the absence of dis- 
sipation the characteristic polynomial of this system is 
Pa = A and the graph has a zero eigenvalue of multiplic- 
ity 4. However, if a small amount of dissipation e is, for 
illustrative purposes, present on vertex 'A' then the char- 
acteristic polynomial becomes p& = /i 4 +e/i 3 —e which has 
roots n k = e 1 / 4 e 2mfe / 4 - e/4 + <3(e 7 / 4 ) for k = 1 ... 4 as 
e — > and therefore /^ max > (to see this set \j! = e -1 / 4 fi 
to give the rescaled problem // 4 + <5// 3 — 1 = 0, where 
S = e 3 / 4 , use the ansatz fjf = /1q+<5// 1 +. . . and solve in the 
limit 5 — > 0) . In this case changes in cyclic structure dur- 
ing the evolutionary process may not necessarily trigger 
concordant changes in stability. For instance removing 
either of the two edges connected to vertex 'A' in Fig. [5] 
decreases $ and "J/ yet increases /Lt max (see main text for 
definitions of $, ^ and /i max )- In practice, since they 
rely on rather specific structural configurations, changes 
in stability which are discordant with changes in cyclic 
structure are rare in our model. Nevertheless, this obser- 
vation highlights the general case: stability is determined 
by the interplay between the cycles of S(t), their sign, 
and the relative rates of decay. 
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